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Abstract. 

Finite difference schemes are here solved by means of a linear matrix equation. The theoretical 
study of the related algebraic system is exposed, and enables us to minimize the error due to a finite 
difference approximation. 

Key words. Finite difference schemes, Sylvester equation 

AMS subject classifications. 15A15, 15A09, 15A23 

1. Introduction: Scheme classes. Finite difference schemes used to approxi- 
mate linear differential equations are usually solved by means of a recursive calculus. 
We here propose a completely different approach, which uses an equivalent matrix 
equation. 



Consider the transport equation: 
fjii fill 

(1.1) ^ + c|^ = , X e [0,L], t e [0,T] 

with the initial condition u{x, t = 0) = uq{x). 
Proposition 1.1. 

A finite difference scheme for this equation can be written under the form: 

n+l , o n I n— 1 . r n , n , /■ n+1 , n — 1 . /i n+1 . n n— 1 n 

aui ^ +pu^ +7Ui +dui+i +eui-i +Cui+i +r]Ui^i +Oui-i ^ +§Ui+i =0 

(1.2) 

where: 

(1.3) ur = u{lh,mT) 

I Cz {i — 1, z, i + 1}, TO G {n — I, n, n + 1}, j — 0, Ux, n = 0, nt, h, T denoting 

respectively the mesh size and time step (L — n^h, T = nt t). 

The Courant-Friedrichs-Lewy number (of I) is defined as a — ct/H . 



A numerical scheme is specified by selecting appropriate values of the coefficients 
a, 13, 7, 8, £, Q, rj, 6 and i? in equation Values corresponding to numerical 

schemes retained for the present works are given in Table \]~7\ 

The number of time steps will be denoted nt, the number of space steps, nx- In 
general, nx^ nt- 
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Table 1.1 
Numerical scheme coefficient. 



Name 
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/3 
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Crank-Nicolson 
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— c 

h7 


— c 





— c 


— c 

717 






The paper is organized as follows. The equivalent matrix equation is exposed in 
section El Specific properties of the involved matrices are exposed in section |2| A 
way of minimizing the error due to the finite difference approximation is presented in 
section 0] 

2. The Sylvester equation. 

2.1. Matricial form of the finite differences problem. Let us introduce the 
rectangular matrix defined by: 



(2.1) 



U — [Ui"] l<i<„^-l, l<n<nt 



Theorem 2.1. The problem J^l.^) can he written under the following matricial 



form: 
(2.2) 



MiU + U M2 + C{U) ^ Mo 



where Mi and M2 are square matrices respectively Ux — ^ by Ux — 1, nt by nt, given 
by: 



(2.3) Ml 



( 13 5 ... \ 

■■• ■•• ■■• 

: ■■• ■■• (3 6 

V ... e P J 

the matrix Mq being given by: 



M2 



/ 7 
a ■•■ 

■■• ■•■ 



\ 



V 



a J 



Mo 



(2.4) 



V -7<_i-5< 



-T,ul 



■riu„ ■ 












-S ul^ — i) u 



and where C is a linear matricial operator which can be written as: 
(2.5) £ = Ci+C2 + C3 + C4 
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where L\, C2, and £4 are given by: 



( 



A(t/)-C 

(2.6) 



V 



\ 



/ 



^2(C/)=r7 




















u\ 











u\ 


u\ 




U2 



V 



(2.7) 



/ 



\ 




J 



Ci{U) = d 



( 




ui 



V 



Proposition 2.2. The second member matrix Mq bears the initial conditions, 
given for the specific value n = 0, which correspond to the initialization process when 
computing loops, and the boundary conditions, given for the specific values i — 0, 
i = Ux- 

Denote by Uexact the exact solution of Hl.l|l . 
The corresponding matrix Uexact wiU be: 



(2.8) Uexact — \Uexacti^]l<i<n^-l,l<n<nt 

where: 



(2.9) Uexacti — Uexacti^'^i^t^i) 

with Xi — i h, tn — nr. 

Definition 2.3. We will call error matrix the matrix defined by: 

(2.10) E = U -Uexact 



Consider the matrix F defined by: 

(2.11) F = Ml Uexact + Uexact M2 + C{Uexact) " Mq 

Proposition 2.4. The error matrix E satisfies: 



(2.12) 



MiE + EM2+C{E) = F 
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2.2. The matrix equation. 

2.2.1. Theoretical formulation. Theorem 2.5. Minimizing the error due to 
the approximation induced by the numerical scheme is equivalent to minimizing the 
norm of the matrices E satisfying ^2.1I1\} . 

Note: Since the linear matricial operator C appears only in the Crank-Nicolson 
scheme, we will restrain our study to the case C — 0. The generalization to the case 
£ ^ can be easily deduced. 

Proposition 2.6. 
The problem is then the determination of the minimum norm solution of: 



(2.13) MiE + EM2^F 

which is a specific form of the Sylvester equation: 



(2.14) 



AX + XB = C 



where A and B are respectively m by m and n by n matrices, C and X , m by n 
matrices. 



The solving of the Sylvester equation is generally based on Schur decomposition: for 
a given square n by n matrix A, n being an even number of the form n = 2p, there 
exists a unitary matrix U and a upper triangular block matrix T such that: 



(2.15) 



A = U*TU 



where U* denotes the (complex) conjugate matrix of the transposed matrix ^U. The 
diagonal blocks of the matrix T correspond to the complex eigenvalues A; of A: 



( ^1 




(2.16) 



T = 







\ 



V 

where the block matrices Ti, i = 1, 







p are given by: 



(2.17) 



TZe [Xi] 2m [Xi] 
— 2m[Xi\ TZe[Xi] 



TZe being the real part of a complex number, and Xm the imaginary one. 

Due to this decomposition, the Sylvester equation require, to be solved, that the 

dimensions of the matrices be even numbers. We will therefore, in the following, 
restrain our study to — 1 and nt being even numbers. So far, it is interesting to 
note that the Schur decomposition being more stable for higher order matrices, it 
perfectly fits finite difference problems. 
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Complete parametric solutions of the generalized Sylvester equation H2.13|l is given in 

12], m 



As for the determination of the solution of the Sylvester equation, it is a major topic 
in control theory, and has been the subject of numerous works (see 0, 0, |H], [H], 

m, m)- 

In the method is based on the reduction of the he observable pair {A,C) to 
an observer-Hessenberg pair {H,D), H being a block upper Hessenberg matrix. The 
reduction to the observer-Hessenberg form (iJ, D) is achieved by means of the staircase 
algorithm (see 0, ...). 

In [^, in the specific case of B being a companion form matrix, the authors propose 
a very neat general complete parametric solution, which is expressed in terms of 
the controllability of the matrix pair {A,B), a symmetric matrix operator, and a 
parametric matrix in the Hankel form. 

We recall that a companion form, or Frobenius matrix is one of the following kind: 



(2.18) 



B 



( 
1 





V 



-6o \ 

-hi 



1 / 



These results can be generalized through matrix block decomposition to a block com- 
panion form matrix: 



(2.19) 



Af2 = 



V 











Mf 




















\ 



the M^^ , 1 < P < being companion form matrices. 



Another method is presented in T4', where the determination of the minimum- norm 

solution of a Sylvester equation is specifically developed. 

The accuracy and computational stability of the solutions is examined in (5|. 



2.2.2. Existence condition of the solution . Theorem 2.7. In the specific 
cases of the Lax and Lax- Wendroff schemes, {2.14\j has a unique solution. 

Proof. H2.14|l has a unique solution if and only if A and B have no common 
eigenvalues. 

The characteristic polynomials Pmi , F'A-h of and M2 , can be respectively calculated 
as the determinants of respectively ^''2' by , ^ by ^ diagonal block matrices: 



(2.20) Pmi(A) = ((A - (3f - 5e)^ , Pm.W = - a^)"^ 
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In the specific cases of tlie Lax and Lax-Wendroff scliemes: 07 = 67^ ±/3 + y5e. 
Hence, (|2.12ll lias a unique solution, which accounts for the consistency of the given 
problem. 

□ 

Corollary 2.8. In the specific case of the Leapfrog scheme, \2.14^ has a unique 
solution if and only if t ^ ^. 

3. Specific properties of the matrices Mi and M2- 

3.0.3. Inversibility of the matrix Mi . Theorem 3.1. In the specific case 
of the Lax and Leapfrog schemes, Mi is inversible. 

Theorem 3.2. In the specific case of the Lax-Wendroff scheme. Mi is inversible 
if and only if 



(3.1) 



-1 , c'r 2^ {a'-iy 



Proof. The determinant Di of Mi can be calculated as the determinant of a 

2~ 



by ^\ ^ diagonal block matrix 



(3.2) 
□ 



Di = iP^-Se)- 



3.0.4. Nilpotent components of the matrix M2. Proposition 3.3. M2 
can be written as: 

(3.3) M2^a'^N + jN 

where N is the nilpotent matrix; N denotes the corresponding transposed matrix: 

/ 1 ... \ 
■■• ■•• : 
: : ■■. ■■. 



(3.4) 



N 



V 



1 
... 0/ 



Proposition 3.4. 
For either a = 0, or 7 = 0, M2 will thus be nilpotent, of order Uf. 

Proposition 3.5. 

For J — 0, if Ml is inversible, the solution at t = nt dt can then be immediately 
determined. 

Proof. 

In such a case, multiplying (|2.2|) on the right side by 1/2"'^^ leads to: 
(3.5) Ml U M2"*"^ + U M2"* = Mo Ma"*"^ 
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1. e.: 

(3.6) 

which leads to: 
(3.7) 



Mii/Ma"'"^ = M0M2 



f/Ma"*"^ = Mf ^ Mo Ms"*"^ 



Due to: 



(3.8 



we have: 



(3.9) 



/ ... \ 

: ■•• ■■• : 

: : ■•• ■•. 

: ■•• 

V 1 ... / 



/ Win, ... \ 

U2nt '•■ ■■■ : 
: : ■•. ■•. 




0/ 



The solution a.t t = rit dt can then be immediately determined. 



(4.1) 



4. Minimization of the error. 
4.1. Theory. Calculation yields: 



/3 (5 + e) + 13^ 



) ( 



fi{5 + e) 
l3{& + e) £2 _^ /32 



-»((r :o (r :o 



)) 



The singular values of Mi are the singular values of the block matrix ( 



P (S + e) £2 + /32 



1. e. 



(4.2) 



of order ^^^'^j and 



(4.3) 



of order 



1 (2/j2 + J2 ^ ^2 ^ ^ ^4^2 + 52+£2_2j£) 



2 ■ 



The singular values of M2 are a^, of order and 7^, of order 
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Consider the singular value decomposition of the matrices Mi and M2: 



(4.4) U^MiVi = j , U^MiV, = ( 

where Ui, Vi, U2, V2, are orthogonal matrices. Mi, M2 are diagonal matrices, the 
diagonal terms of which are respectively the nonzero eigenvalues of the symmetric 
matrices Mi^Mi, M2^M2. 

Multiplying respectivelv 12 . 1 3l on the left side by '^Ui, on the right side by V2, yields: 



(4.5) t/f Ml EV2 + Ul E M2 V2 = Ul F V2 
which can also be taken as: 

(4.6) '^Ui Ml Vi^ViE V2 +^ J7i S ^f72 ^J72 M2 V2 = F V2 
Set: 



(4.7) ^ViEV2^l ^]jUiE^U2^l §3 

V ^21 E22 J \ E21 E22 



(4.8) 'U1FV2 
We have thus: 



Ttt tpm ^ I ^ ^ 
F21 F22 



(4.9) 
It yields: 



Ml Ell M1E12 \ , E11M2 \ / Fii F] 



12 



^ ^ ' \ E21M2 Q j \ -^21 F22 



Ml Eii+Jhi M2 = Fn 
(4.10) { NhEi2 = F12 

E21 M2 = F21 

One easily deduces: 
(4.11) 

The problem is then the determination of the En and £^11 satisfying: 



E12 = Ml ^ F12 
E21 = F21 M2 



(4.12) 
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Denote respectively hy efj, efj the components of the matrices E, E. 
The problem 14 . 1 21 uncouples into the independent problems: 
minimize 



(4.13) 



-2 ^2 



under the constraint 

(4.14) Mu.^Tj + Ah^^^j ^fTuj 

This latter problem has the solution: 



(4.15) 



T i 



The minimum norm solution of 12.131 will then be obtained when the norm of the 
matrix Fn is minimum. 

In the following, the euclidean norm will be considered. 
Due to (gSJl: 

(4.16]|^^|| < < \\Ui\\ \\F\\ \\V2\\ < llC/lll IIF2II \\MiUe^act + Ue^actM2-Mo\\ 

Ui and V2 being orthogonal matrices, respectively n^; — 1 by n^; — 1, by nt, we have: 

(4.17) ||[/i||2^n,-l , \\V2f^nt 
Also: 

(4.18) ||Afi||2 = !^(2/32 + j2_^e2) , \\M,f = !^ + j^) 

The norm of Mq is obtained thanks to relation (|2.4() . 
This results in: 

ll^^ll < Vnt (n, - 1) |||[/e.act|| {\I^^VW+S^T1^+^V^^^T^) + IIM 
(4T9) 

||i^ii|| can be minimized through the minimization of the right-side member of H4.19|l . 
which is function of the scheme parameters. 

4.2. Numerical example: the specific case of the Lax scheme. For the 

Lax scheme: 7 = 0. The remaining coefBcients (see Table II. Ill can be normalized 
through the following change of variables: 

a = ha 

e = he 



'01 
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Set: 



(4.21) 
Thus: 



Afo ^ hMo 



(4.22) ||Mo||2 = r5^. 

n=l 

Advect a sinusoidal signal 



nt 

2 , ^2 '-2 



(4.23) 



. 2 vr 

u = cos [ — — [X — ct) \ 
A 



through the Lax scheme, with Dirichlet boundary conditions. 
The right-side member of H4.19|l is then: 



1 1 



' 2 2 ^1 1 



2 2cfl 



nt ul + 



n.r, - 1 



1 1 



1 1 



2 2cflJ \2 2cflJ \/2 



2 2cfl 
(4.24) 

where uq, uj^ respectively denote the Dirichlet boundary values at x = 0, x — L. 
It is minimal for the admissible value cfl = 1. 

The value of the L2 norm of the error, for two significant values of the cfl number 
(Case 1: cfl = 0.9; Case 2: cfl — 0.7), is displayed in Figure Wa\ The error curve 
corresponding to the value of the cfl number closest to 1 is the minimal one. 



1 . 2 

1 

0.8 
0.6 
0.4 
0.2 



20 



40 60 80 100" 



n 



Case 1 
Case 2 



Fig. 4.1. Value of the L2 norm of the error for different values of the cfl number. 



5. Conclusion. Thanks to the above results, we here propose: 

1 . to study the intrinsic properties of various schemes through those of the re- 
lated matrices Mi and M2; 

2. to optimize finite difference problems through minimization of the symbolic 
expression of the error as a function of the scheme parameters. 
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